#include <string.h>

#include "lapack.h"
#include "BandedSystem.hpp"

BandedSystem::BandedSystem(int n, int kl, int ku) :
  m_n(n), m_kl(kl), m_ku(ku), m_ldab(2*kl+ku+1),
  m_ldb(n), m_info(0) {
  m_AB = new double[n*m_ldab];
  m_B = new double[n];
  m_ipiv = new int[n];
	memset(m_AB, 0, n*m_ldab);
	memset(m_B, 0, n);
	memset(m_ipiv, 0, n);
}

BandedSystem::~BandedSystem() {
  delete[] m_AB;
  delete[] m_B;
  delete[] m_ipiv;
}

int BandedSystem::LUDecompose() {
  dgbtrf_(&m_n, &m_n, &m_kl, &m_ku, m_AB, &m_ldab, m_ipiv, &m_info);
  return m_info;
}

int BandedSystem::solveLU() {
  char trans = 'N';
  int nrhs = 1;
  dgbtrs_(&trans, &m_n, &m_kl, &m_ku, &nrhs, m_AB, &m_ldab, m_ipiv,
          m_B, &m_ldb, &m_info);
  return m_info;
}
